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THE STOCHASTIC EVOLUTION OF A PHOTOCELL. 
THE GILLESPIE ALGORITHM IN A DYNAMICALLY VARYING VOLUME. 



T. CARLETTI AND A. FILISETTI 



Abstract. In the present paper we propose an improvement of the Gillespie algorithm allowing 
us to study the time evolution of an ensemble of chemical reactions occurring in a varying 
volume, whose growth is directly related to the amount of some specific molecules, belonging to 
the reactions set. 

This allows us to study the stochastic evolution of a protocell, whose volume increases because 
of the production of container molecules. Several protocells models arc considered and compared 
with the deterministic models. 

1. Introduction 

All known life forms are composed of basic units called cells; this holds true from the single-cell 
prokaryotc bacterium to the highly sophisticated eucaryotes, whose existence is the result of the 
coordination, in term of self-organization and emergence, of the behavior of each single basic unit. 

While present day cells are endowed with highly sophisticated regulatory mechanisms, that 
represent the outcome of almost four billion-years of evolution, it is generally believed that the 
first life-forms were much simpler. Such primordial life-bricks, the protocells, were most probably 
exhibiting only few simplified functionalities, that required a primitive embodiment structure, a 
protometabolism and a rudimentary genetics, so to guarantee that offsprings were "similar" to 
their parents mfTSlfTT). 

Intense research programs are being established aiming at obtaining protocells capable of growth 
and duplication, endowed with some limited form of genetics }12[ [T51 [HI [T7j . Despite all efforts, 
artificial protocells have not yet been reproduced in laboratory and it is thus extremely impor- 
tant to develop reference models |31 UHl [111 US] that capture the essence of the first protocells 
appeared on Earth and enable to monitor their subsequent evolution. Due to the uncertainties 
about the details, high-level abstract models are particularly relevant. Quoting Kaneko [7] it is 
necessary to consider "simplified models able to capture universal behaviors, without carefully 
adding complicating details" . 

Most of the models present in the literature are based on deterministic differential equations 
governing the evolution of the concentrations of the involved reacting molecules. Even if the 
results are worth discussing and provide important insights, it should be stressed that the former 
assumptions are rarely satisfied in a cell [S]. Firstly, the number of involved molecules is small 
and should be counted by integer numbers, hence the use of the concentrations can be questioned; 
secondly, the presence of the thermal noise introduces in the system a degree of stochasticity than 
cannot be trivially encoded by a differential equation, mostly because this makes the time evolution 
a stochastic process. One possible way to overcome such difficulties is to use the Chemical Master 
equation: given the present state of the system, namely the number of available molecules for each 
species, and the possible reactions among them, one can compute the transition probabilities to 
reach and leave the given state and thus get a partial differential equation describing the time 
evolution of the probability distribution of having a given number of molecules at any future 
times [3 [6]. 

Analytically solving the resulting equation is normally a very hard task, one should thus resort 
to use numerical methods. A particularly suitable one is the algorithm presented by Gillespie 0(5], 
allowing to determine, as a function of the present state of the system, the most probable reaction 
and the most probable reaction time, i.e. the time at which such reaction will occur. 
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Let us however observe that in the setting we are hereby interested in, the chemical reactions 
occur in a varying volume, because of the protocell growth; we thus need to adapt the Gillespie 
method to account for this factor. To the best of our knowledge, there arc in the literature very 
few papers dealing with the Gillespie algorithm in a varying volume [H |9] . Moreover in all these 
papers, the volume variation can be considered as an exogenous factor, not being directly related to 
the number of lipids forming the protocell membrane. So our main contribution is to improve the 
Gillespie algorithm taking into account the protocell varying volume which is moreover consistent 
with the increase of the number of lipids constituting the protocell membrane. 

The paper is organized as follow. In Section [2] we briefly recall the Surface Reaction Models 
of protocell, that would be used to compare our stochastic numerical scheme. Then in Sections [3] 
and [4] we will present our implementation of the Gillespie algorithm in a dynamically varying 
volume. Finally in Section [5] we will present some applications of our method. 

2. Surface Reaction Models 

Among the available models for protocells, a particularly interesting one is the Surface Reactions 
Model [21 [in], SRM for short, and its applicability to the synchronization problem. Such model 
is roughly inspired by the Los Alamos bug hypothesis [131 [U] but which, due to its abstraction 
level, the SRM can be applied to a wider set of protocell hypotheses. 

The SRM is build on the assumption that a protocell should comprise at least one kind of 
"container" molecule (typically a lipid or amphiphile), hereby called C molecule, and one kind 
of replicator molecule - loosely speaking "genetic material", hereafter called Genetic Memory 
Molecule, GMM for short, and named with the letter X . There are therefore two kinds of reactions 
which are crucial for the working of the protocell: those which synthesize the container molecules 
Eq. ([IJ and those which synthesize the GMM replicators Eq. ([2]) 

(1) X, + L,^^X, + C, 



(2) X, + P, >X, + X, . 

In both cases Li and Pj are the buffered precursors, respectively of container molecules and of the 
j"th GMM, while Ui and Mij are the reaction kinetic constants. 

A second main assumption of the SRM, is that such reactions occur on the surface of the 
protocell, exposed to the external medium where precursors are free to move. Hence, as long as 
container molecules are produced, they are incorporated in the membrane that thus increases its 
size, until a critical point at which, due to physical instabilities, the membrane splits and two 
offsprings arc obtained, each one getting half of the mother's GMMs and whose size is roughly 
half that of the mother just before the division. 

Under the previous assumptions and in the deterministic setting, one can prove [3l I16j that the 
number of membrane molecules and the number of GMMs evolve in time according to: 

\c -(^r<^-x 

I X = C^-^M ■ X , 



(3) 



where X = {Xi, . . . ,Xn) represents the amount of each GMM, a = (ai,...,aAr) is the vector 
of the reaction constants responsible for the production of C molecules from the X molecules 
plus some appropriate precursor. (Mij) denotes the reaction constant at which Xi is produced 
by Xj plus some precursor. (3 € [2/3,1] is a geometrical shape factor that relates the surface to 
the volume of the protocell and p is the lipid density (for more details the interested reader can 
consult (SI [Ml). Let us observe that in this setting the precursors are assumed to be buffered and 
thus their amount to be constant, hence the latter can be incorporated into the constants a and 
M. 

So starting with a initial value of container molecules, C{to) = Co, and of GMMs, X{tQ) = Xq, 
the protocell will grow until some time to + ATi at which the amount of C molecules has doubled 
with respect to the initial value, C(to + ATi) = 2Co and thus the protocell undergoes a division. 
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Each offspring will get half of the GMMs the mother protocell had just before the division, X^^^ = 
X{to + ATi)/2. And the protocell cycle starts once again. One can prove [3l[T6| that under suitable 
conditions X*-"^ tends to a constant value once n goes to infinity, implying thus the emergence of 
synchronization of growth and information production. 

3. The method 

Let us now improve the previous scheme by introducing a probabilistic setting a la Gillespie. 
We thus consider a protocell made by a lipidic vesicle and containing a well stirred mixture of N 
GMMs, Xi, . . . ,Xn, that may react through m elementary reaction channels i?^, fi = 1, . . . ,to, 
running within the volume V{t) of the protocell. 

Let us observe that because of the protocell growth the volume is an increasing function of 
time. Actually one can relate the volume to the amount of container molecules via their density 
V = C/p where C denotes the integer number of molecules forming the lipidic membrane. We 
will hereby use the same symbol Xi to denote both the i-th GMM and the integer number of 
molecules of type Xi in the system. 

For each reaction channel i?^ assume that there exists a scalar rate such that c^^dt + o{dt) is 
the probability that a random combination of molecules from channel i?^ will react in the interval 
[t,t + dt) within the volume V{t). 

Let h^{Y) be the total number of possible distinct combinations of molecules for a channel 
when the system is in state Y = {Xi, . . . , Xn,C), then we can define the propensity ^ of the 
reaction i?^ to be a^{Y) = /i^(Y')c^. 

One can prove [5] that for a binary reaction the rate can be written in the form = kf^jV , 
where is a fixed constant. Similarly one can prove that for a reaction involving n different 
species, we get: = And thus for a single molecule reaction, i.e. a decay, we get 

= fc^, namely independently from the volume. 

Let us now assume that among the m reactions, Qi involve one single molecule, Q2 are binary 

reactions, Q3 are ternary reactions and so on. Of course Qi + (52 ^ ^ Qn+i = rn. We recall that 

we have N GMMs and the container type molecule C, hence + 1 species. For short we will 
denote Qi the set of indices jj, for mono molecule reactions, and by Q the remaining ones. Let us 
observe that in this way some coefficient a^, will depend both on the system state Y and on the 
time via the volume V{t): af^{Y,t) for p e Q. 

More precisely to study the time evolution of the system we need to determine the probability 
Pfj,{T\Y,t)dT, that given the system in the state Y = (Xi, . . . , X„, C) at time t, then the next 
reaction will occur in the infinitesimal time interval (t + T,t + t + dr) and it will be the reaction 
R^. This probability will be computed as 

(4) P^{T\Y,t)dT = Pr,ot{T\Y,t)xaf,{Y,t + T)dT, 

where Pnot{T\Y,t) is the probability that no reaction occurs in (t,t + T) given the state Y at time 
t, whereas the rightmost term denotes the probability to have a reaction R^ in (t + T,t + t + dr) 
given the state Y at time t + t. 

To compute the first term P„ot, let us take s e + r] and observe that: 

Pnotis + ds\Y, t) = Pnot{s\Y, t)Pnot{ds\Y, t + s) = Pnot{s\Y, t) ^1 - a^,{Y,t + .s)ds j , 

being 1 - Z^ct/^C^, ^ + s)ds the probability that no reaction will occur in (t + s,t + s + ds) once 
we are in state Y at time t + s. Thus rewriting the previous difference equation as a differential 
equation, passing to the limit ds 0, and observing that Pnot{0\Y,t) = 1, we get the solution: 



(5) P„ot(T|r,i) =exp 

where 



-Aq,{Y)t- JJ AQ{Y,s + t) ds 



AQdy)=Za^(Y) and Aq{Y,s + t) = a^{Y,s + t) . 
The apparent asymmetry in the exponential term in ([5]) is easily recovered by observing that 

AqAy)t = i;; AQ,{Y)ds. 
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We can thus conclude that 



(6) P,,(T|r,t)dr = exp 



Xt+T 



ds 



a^,{Y,t + T)dT . 



Let us observe that the rightmost term is correctly a^j,{Y,t + r), namely the system is still in the 
state Y at time t + t, because no reaction has been produced in (t, t + t). 

Let us recall that the volume enters in the previous relation via the function Aq , more explicitely 
one has 



(7) 



Aq{Y,s)= 



h^{Y)k, 
V{s) 



E 



h^{Y)k^ 
{V{s)Y 



E 



K{Y)k^ 

{v{s)r 



that can be rewritten in terms of C molecules using the relation C = pV . So our method applies 
to a different problem with respect to the one considered in [5], in fact in our case the volume 
growth is not imposed a priori but dynamically evolves accordingly to the reaction scheme, if C is 
produced then V increases otherwise it will keep a constant value, while in [3] the volume growth 
is an exogenous variable. 



4. The stochastic simulation algorithm in a growing volume 

Once we have the probability function P^{T\Y^t) we can build an algorithm that reproduces 
the time evolution given by the model defined above. 

Given the system in some state Y at time t, we must determine the interval of time t and 
the reaction channel i?^ according to the probability distribution function P^{T\Y,t), and finally 
update the state Y ^ Y + v^, where is a stoichiometric vector representing the increase and 
decrease of molecular abundance due to the reaction R^^. This will be accomplished following the 
standard approach by Gillespie [5] but taking care of the time dependence of the propensities. We 
will thus need to compute the cumulative probability distribution function and then make use of 
the inversion method [5], to determine the channel ^ and the next reaction time t, distributed 
according to P^,{T\Y,t). 

From we can compute the cumulative distribution function 



(8) 



Jo n 



providing the probability that any reaction will occur in (t,t + T) starting from the state Y at time 
t. The function F{t\Y, t) can be explicitely computed by 



Proposition 4.1. Under the above assumptions we have 



(9) 



F{T\Y,t) = 1-exp 



-Aq,{Y)t- f£ Aq{Y,s) 



ds 



Proof. The first step is to use ^ and perform a sum over all the channels iJ. to rewrite (jS]) as 



F{T\Y,t) = J\AQ,{Y)+AQ{Y,t + s))cxp]^-AQ,{Y)s- AQ{Y,r)dr 



ds . 



Then we can observe that 
d 



and thus 



— (cxp[-AQ,(y)s-^ "AQ(r,r)dr]) = 
= - {Aq, (Y) + Aq{Y, t + s)) exp [-Aq, {Y)s- Aq{Y, r) dr 

F{T\Y,t) = - JJ j-(^exp^-AQ,{Y)s- AQ{Y,r)drJjds 
= 1-cxp^-Aq,{Y)t- j£ AQ{Y,r)dr 



□ 
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Once we have the cumulative distribution function we can obtain the value r by drawing a 
radom number ui from an uniform distribution in [0, 1] and then solve with respect to r the 
implicit equation: 

(10) ui = 1- exp 



-Aq,{Y)t-J^ ^ AQ{Y,s)ds 



Let us stress once again that this is not as straightforward as for the original Gillespie [5] scheme, 
or the simplified one presented in [9] , because of the time dependence of Aq via the volume. One 
can nevertheless found suitable approximation for the integral, this will be the goal of the next 
sections. 

4.1. The adiabatic assumption. Let assume that r is very small, or which is equivalent, that 
the time scale of the chemical reactions involving the GMMs is much faster than the production 
of container molecules, hence the volume growth is very slow compared with the production of the 
chemicals Xi. 

Under this hypothesis one can assume that in the interval {t,t + T) the volume doesn't vary and 
thus one can made the following approximation 

(11) AQ{Y,s)ds^AQ{Y,t)T. 

One can thus explicitely solve equation (fTO]) to get: 



that is the standard Gillespie result except now that AQ{Y,t) depends on time and as long the 
volume increases, then the contribution arising form AQ{Y,t) mights become smaller because 

Aq ~ 1/V. 

4.2. The next order correction. One can obtain a somehow better estimate valid in the case 
of comparable time scales for the reactions involving GMM and the container growth. The idea 
is to compute the integral in Eq. ()10|) using the following approximation: 



f^''^ AQ{Y,s)ds = JJ AQ{Y,t + s)ds= fj(^AQ{Y,t) + ^^^^s + ...^ds 
(13) - AQiY,t)rJ-M^I^,oir^). 

where dAQ{Y,t)/dt can be obtained using the definition ([7]) and expressing the volume in terms 
of C = V{t)p, namely: 

dAQ{Y,t) _ C'/ ^ Mm^^2 ^ Mmi^...^^ ^ h,{Y)k^\ 



To compute C/C we make the assumption that in a very short time interval, as the one we are 
interested in, the deterministic growth of the container is a good approximation for the stochastic 
underlying mechanism; this implies that we can use ([3]) 

C (C{t)\^'^ a-X{t) 



\ P ) i 



C \ p I C{t) 

Inserting the previous result into (|13l) and finally solving (|10p with respect to r, we can compute 
the next reaction time up to correction of the order of r^, as: 



-(Aq,(Y) * AQ(Y,t)) + J{Aq,(Y) + AQ(Y,t)y -2\og(l -«i)iq(r,f) 

A^) • 

where we wrote for short Aq{Y, t) = dAQ{Y, t)/dt and we selected the positive square root in such 
a way in the limit AQ{Y,t) ^ we recover the previous solution (|12p . 
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Remark 4.2 (On the existence of tgui)- In the case of variable volume a new phenomenon can 
arise: the volume growth can he so fast that no reaction can occur in the interval (t,t + T + dr) for 
any t. Mathematically this translates into a sign condition for the term under square root in 
if-- 

(15) log(l-ui) < {AQ,{Y)+AQ{Y,t)f/{2AQ(Y,t)), 

then equation (jlOp has no real solution. 

This can be geometrically interpreted as follows. The relation (jlOp determines tqui as the 
intersection of the parabola -Aq-^{Y) - AqiY^t^r - AQ{Y,t)T^ /2 with the horizontal line log(l - 
ui), which is negative because ui € (0,1). Such parabola intersect the y-axis at ti = and T2 = 
-2{Aq-^(Y) + AqiY^t)) I AQ{Y^t) > and it is concave. Then its absolute (negative) minimum 
is reached at the vertex Ty = {ti + t2)/2 and its value is {Aq-^{Y) + Ag(Y, t))^/(2AQ(Y, t)) and 
it is negative because AQ(Y,t) is negative. Hence if the horizontal line is below this value, i.e. 
condition (jlSp is verified, the parabola and the line do not have any real intersections (see Fig. Qp. 




Figure 1. Geometrical interpretation of the existence of the next reaction time 
Tgui. Left panel : tgui is the smallest intersection between the parabola and the 
horizontal line log(l - ui). Right panel tgui doesn't exist, the horizontal line is 
located below the minimum of the parabola. 

Let us also observe that, whenever it exists, tgui is always positive as it should be. In the case 
of a protocell the non existence of such next reaction time could be translated into the death by 
dilution of the protocell. 

4.3. The next reaction channel. Whenever the next reaction time does exist, the next reaction 
channel is determined using the classical Gillespie method, namely by drawing a second uniformly 
distributed random number U2 g [0, 1] and fix the channel fi such that: 

(16) Z''-(Y,t + T)<U2ao{Y,t + T)< Y,a,{Y,t + T), 

where ao{Y, t + T) = Aq, {Y) +AQ{Y,t + T) = YZi a^Y. t + r). 

Remark 4.3. Let us observe that if all the reactions involve the same number of chemicals, 
then the determination of which reaction channel fi will be activated in the next reaction, doesn't 
depend on the volume which factorizes out from (jl6p . In fact assuming all the reactions to involve 
p chemical, we obtain by definition 

a^{Y,t + T) = —- -p Vz.6{l,...,m}, 

[V{t + T)Y 

and thus (fT6l) rewrites: 



't^ h^{Y)k^ " K{Y) 

h mt + T)r - h mt + T)r - h mt + T)r ' 

which is clearly independent of the volume value V. 
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5. Some applications 

The aim of this section is to provide some appHcations of the previous algorithm to the study 
of the evolution of a protocell. 

5.1. One single Genetic Memory Molecule. The simplest model is the one where only one 
GMM specie is present in the protocell [16] and thus only two chemical channels are active: 

■n 

channel 1, i?i : X + Pi >2X 

a 

(17) channel 2, i?2 : X + >X + C, 

where Pi and Li arc, respectively, precursors of GMM, i.e. nucleotide, and precursors of am- 
phiphiles. 

One can thus compute the propensities in the state Y = {X, C) at time t: 

(18) ai{X,C,t) = hi(X,C)-^ = r]^-^ and a2(X,C,t) = h2(X,C)-^ = ; 

let us observe that we assume that precursors arc buffered and thus they are constant. 

Because system (fTT)) contains only bimolecular reactions, all the propensities are time depen- 
dent, hence Aq^ = and Aq = ai{X,C,t) +a2{X,C,t) = {Pirj + Lia)X /V{t), thus (UHl) simplifies 
into ^ 

ui = l-cxp|^-^ AQ{Y,s)ds , 
whose second order solution (|14p is given by 



-AQ{Y,t)+ ^J{AQ{Y,t))^ -2\og{l -ui)AQ{Y,t) 



AQ{Y,t) 
and 

dAqiX, C,t) _ V{t) j PiyX LiaX\\ IC\^^^ pLiaX'^ 

dt " V{i^\v{t) ^ V{t) )\v{t)=c{t)ip~ \7/ C2 

So we can finally obtain 



{Pirj + Lio) 



^""'^ LiaXXc) \[LiaX\c) 
provided 

log(l-ui)>-^(^)'"' ^'{Piv + Lia). 
Which reaction channel yu will active in the time interval [t, t + r] can be obtained according to 

. {PiV + Lia)X ^ Pir^X , p,^ . 

ifu2 < — ^ namely < ^2 < p^,,+^^ then /i = 1 

..PivX {Piri + Lia)X (Piri + Lia)X , p,„ 

if^-!—<U2^^^ LJ_<.L_L! LJ— namely ^ < M2 < 1 then u = 2 . 

Let us observe that according to remark 14.31 the choice of ^ doesn't depend on the volume, 
because only binary reactions are present. 

Let Co be the initial amount of container molecules, then we assume that once C{t) = 2Cq the 
protocell splits into two offspring, almost halving the GMM amount. More precisely we assume 
that the first offspring will get a number of GMMs drawn according to a Binomial distribution 
with parameter p = 1/2 and n = X{t). From this step, for technical reason, only one randomly 
chosen offspring will be studied during each generation. 

In Fig. [2] we report a comparison between the deterministic ([3|) and the stochastic dynamics, 
under the adiabatic assumption for tghu corresponding to the continuous growth phase of the 
container between two successive divisions. As one should expect, a system composed by a large 



(J2 

+ 2 TT, r log(l - ui) , 

LiapX^{Piri + Lia) 



v2(/3-l) 
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number of molecules exhibits small stochastic fluctuations whose average is not too far from the 
dynamics described by the deterministic model. 




Figure 2. Stochastic vs ODE SRM protocell Case of one GMM, left panel 
the time evolution of the amount of GMM, right panel the time evolution of the 
amount of C. Parameters are : ?/ = 1, a = 1, LI = 500, PI = 600, Xi{0) = 100, 
C(0) = 1000, p = 200 and /3 = 2/3. 



In Fig.[3]we report the amount of GMM, X^'^^ (panel a), at the beginning of each protocell cycle 
and the duplication time (panel b), namely the interval of time needed to double the amount of C 
molecules, for both the stochastic and deterministic models. Once again one can clearly observe the 
small fluctuations of the stochastic system around the value obtained by the numerical integration 
of the deterministic description, Eq. Let us observe that these fluctuations are due to the 
stochastic integrator scheme and also on the division mechanism. 
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Figure 3. Stochastic vs ODE SRM protocell ©. Case of one GMM, left panel 
the amount of GMM at the beginning of each division cycle, right panel the 
division time as a function of the number of elapsed divisions. Parameters are : 
= 1, a = 1, LI = 500, PI = 600, Xi(0) = 100, C(0) = 1000, p = 200 and /? = 2/3. 



We arc now interested in studying the fluctuations dependence on the amount of molecules. 
We already know that for a sufficiently large number of molecules the stochastic dynamics follows 
closely the deterministic one and thus the fluctuations are small. On the other hand, one should 
expect that when the number of molecules decreases, then the fluctuation will rise and the system 
behavior could not be completely described by means of a deterministic approach. This is con- 
firmed by Fig. 2] and Fig. [51 where we can observe that a model composed by a small number of 
initial molecules, 20 times lesser than in the model presented in Fig. [5] exhibits larger stochastic 
fluctuations. 
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Figure 4. Stochastic vs ODE SRM protocell ©. Case of one GMM, left panel 
the time evolution of the amount of GMM, right panel the time evolution of the 
amount of C. Parameters are : rj = 1, a = 1, LI = 500, PI = 600, ^i(O) = 5, 
C(0) = 50, p = 200 and /3 = 2/3. 




20 40 60 SO WO 20 40 60 80 100 



generation number (k) generation number (k) 

Figure 5. Stochastic vs ODE SRM protocell ©. Case of one GMM, left panel 
the amount of GMM at the beginning of each division cycle, right panel the 
division time as a function of the number of elapsed divisions. Parameters are : 
r]=l, a=l, Ll = 500, PI = 600, Xi{0) = 5, C(0) = 50, p = 200 and /3 = 2/3. 



In Fig. |6]we summarize the results of several protocell models each one with a different amount 
of initial molecules, in order to appreciate the influence of the latter on the stochastic fluctuations. 
To compare with, we also report the case of the deterministic model. Because the kinetic constants 
are kept constant, the analytical theory for the deterministic model ensures that the division time 
doesn't vary [3]. Nevertheless the fewer is the initial amount of Xq and Cq, the larger are the 
fluctuations present in the stochastic integration. 



To get a more complete understanding of the fluctuations dependence, we decided to measure 
them using the standard deviation of the protocell division time (after a sufficiently long transient 
phase). In Fig. [7] we report the standard deviation of the division time AT as a function of the 
initial amount of molecules. As expected the fluctuations strength decreases rapidly as soon as 
the number of molecules increases and the relation can be very well approximated by a power law 
distribution with exponent -0.54 ±0.03 (linear best fit). 



5.2. Two non— interacting Genetic Memory Molecules. A slightly more sophisticated model 
can be obtained by considering two linear non interacting GMMs. The system can be described 
by the following chemical reactions: 
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^ X^=50. C„=500 
^X=100. C=1000 
■^ODE 



40 60 80 
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Figure 6. Fluctuation dependence on the initial conditions. We report the divi- 
sion times as a function of the number of elapsed divisions, for 5 different protocells 
models. Protocell o : Xi(0) = 5, C(0) = 10, protocell □: Xi(0) = 10, C(0) = 100, 
protocell v: Xi{0) = 50, C(0) = 500, protocell A: Xi(0) = 100, C(0) = 1000. The 
black line denotes the deterministic protocell. All the remaining parameters have 
been fixed to: r]=l,a = 1, LI = 500, PI = 600, p = 100 and ^ = 1. 
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Figure 7. Fluctuation dependence on the initial conditions. We report the stan- 
dard deviation of the protocell division time as a function of the initial amount of 
molecules Xq (•) and a linear best fit, whose slope is = -0.54 ± 0.03. Parameters 
are: X{0) = 2" with n = 0, 10, C(0) = lOX(O), ?? = 1, a = 1, LI = 500, PI = 600, 
p= 100 and /3= 1. 



channel 1, Pi : Xi + Pi ^2Xi 

ai 

channel 2, P2 : Xi+ Li >Xi + C 

m 

channel 3, P3 : X2 + P2 >2X2 

012 

channel 4, P4 : X2 + P2 >--'^2 + C , 



where P, and Li are, respectively, precursors of the i-th GMM, i.e. nucleotide, and precursors of 
amphiphiles used by the i-th GMM to build a C molecule. 

As previously done, we compare the stochastic and the deterministic models. Results are 
reported in Figure [8] and one can still observe that in presence of a large number of molecules the 
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deterministic dynamics well approximates the stochastic model. On the other hand, the protocell 
division time exhibits large fluctuations around the deterministic value even in presence of a quite 
large number of molecules (sec right panel Fig. IH]). 

The parameters have been set in such a way only one GMM will survive according to the 
analytical theory for the deterministic model. One can observe that, despite the fluctuations, the 
same fate is obtained for the stochastic model (see right panel Fig. [9]). 




Figure 8. Stochastic vs ODE SRM protocell ©. Case of two GMMs, left panel 
the time evolution of the amount of GMM during a division cycle, right panel 
the time evolution of the amount of C molecules. Parameters arc : 771 = 772 = 1, 
ai = aa = 2, LI = 500, L2 = 600, PI = 600, P2 = 670, Xi{0) = ^2(0) = 100, 
C(0) = 1000, p = 200 and /3 = 2/3. 




generaiion number fk) generation number (k) 

Figure 9. Stochastic vs ODE SRM protocell ©. Case of two GMMs, left panel 
the amount of GMM at the beginning of each division cycle, right panel the 
division time as a function of the number of elapsed divisions. Parameters are : 
r]i = ri2 = 1, ai = 02 = 2, LI = 500, L2 = 600, PI = 600, P2 = 670, Xi{0) = ^2(0) = 
100, C(0) = 1000, p = 200 and /? = 2/3. 



Once we reduce the number of involved molecules, the stochastic fluctuations dramatically 
increase (see Fig. [TOl and Fig. [TT|) . 



As in the case of only one GMM, when two non interacting linear GMMs are present the size 
of the stochastic fluctuations as a function of the initial number of molecules follows a power law 
distribution with exponent -0.51 ± 0.05 (linear best fit), sec Fig. [12] the fewer arc the molecules 
in the system, the larger are the fiuctuations around the deterministic dynamics. 
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Figure 10. Stochastic vs ODE SRM protocell ©. Case of two GMMs, left panel 
the time evolution of the amount of GMM during a division cycle, right panel the 
time evolution of the amount of C molecules. Parameters are : ?7i = 772 = 1, 
ai = a2 = 2, LI = 500, L2 = 600, PI = 450, P2 = 670, Xi{0) = XaCO) = 5, 
C(0) = 50, p = 200 and /3= 2/3. 




Figure 11. Stochastic vs ODE SRM protocell ©. Case of two GMMs, left 
panel the amount of GMM at the beginning of each division cycle, right panel 
the division time as a function of the number of elapsed divisions. Parameters 
are : tji = r]2 = 1, ai = a2 = 2, LI = 500, L2 = 600, PI = 450, P2 = 670, 
Xi{0) = X2{0) = 5, C(0) = 50, p= 200 and /3 = 2/3. 



A new phenomenon arises in the case of two GMMs modeled by a stochastic process. There 
can be a breaking of the symmetry emerging in systems composed of two identical GMMs (i.e 
equal kinetic constants, equal initial amounts and availability of precursors) present with a few 
initial amount of each one. Although adopting a deterministic approach the dynamics of the 
two replicators would be perfectly the same, a small fluctuation in the very first instants of the 
protocell evolution entails the dilution of one of the two replicators and thus a different fate for the 
protocell. Let us observe that the probability to have a large fluctuation is never zero, thus waiting 
for a sufficiently long time, a specie can always disappear from the system and thus giving rise 
to the the breaking of the symmetry phenomenon. See Fig. [T3] where we report, as a function of 
the initial amount of molecules Xj(0), i = 1,2, the proportion of simulations where the symmetry 
breaking has been observed repeating 50 times each simulation with the same set of parameters 
and initial conditions during 100 generations. 



6. Conclusion 



In this paper we presented a new stochastic integration algorithm based on the one introduced 
by Gillespie. Our contribution is devoted to the explicit introduction of the volume variation in 
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Figure 12. Fluctuation dependence on the initial conditions. We report the 
standard deviation of the protocell division time as a function of the initial amount 
of molecules Xi(0), i = 1, 2, (•) and a linear best fit, whose slope is = -0.51 ± 0.05. 
Parameters are: Xi(0) = ^2(0) = 2" with n = 0,...,10, C(0) = lOXi(O), 7?i = 772 = 
1, ai = a2 = 2, LI = 500, L2 = 500, PI = 500, P2 = 600, p = 100 and ^ = 1. 
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Figure 13. Symmetry breaking phenomenon. Each point denotes the frac- 
tion of runs exhibiting the symmetry breaking phenomenon, during 100 gen- 
erations, over 50 identical replicas. Parameters are: ^i(O) = ^2(0) = 
[1,2,3,4,5,6,7,8,9,10,25,50], C(0) = lOX, rji = r]2 = 1, ai = 02 = 2, LI = 
500, L2 = 500, PI = 600, P2 = 600, p = 100 and ^ = 1. 



the algorithm, which moreover is directly related to the amount of contained molecules, and thus 
it evolves in a self-consistent way. 

This algorithm straightforwardly adapts to the study of the evolution of a protocell, simplified 
form of cells, where an ensemble of chemical reactions occurs in a varying volume, the volume of 
the protocell, that in turn increases because of the production of container molecules. 

We presented several protocell models and we compare them with the analogous deterministic 
protocell models, namely solved using the ODE. In this preliminary study, we emphasized the 
role of the fluctuations and their dependence on the initial amount of molecules. The dynamics 
is richer than the deterministic one and thus it is worth studying, in particular we deserve to 
future investigations the case where several molecules interact in a linear way but including cross 
catalysis, i.e. the interaction matrix is not diagonal, or they interact in a non-linear way. Also 
the study of the emergence of time-periodic patterns due to the fluctuations, will be analyzed. 
An analytical treatment of the latter case could be possible using some recent technics developed 
by [Hill], see also [2] where the space is also taken into account. 
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